How large is 1ha in size? Below is an example of a 1ha plot at Moore Reef. Note that this is only for demonstration purposes:
Inset is a red box 10 x 10m area that contains 5000 randomly seeded corals in year 1.
This is a smaller subset of a realistic seeding effort: 500,000 individuals across a 100 x 100m 1ha area.
A (very) simple growth calculation (2.5cm per year) is applied to the corals for 5 years with 20% survival
To be updated with coraldynamics
# data packages
library(tidyverse)
library(units)
library(foreach)
# spatial packages
library(sf)
# mapping packages
library(tmap)
library(leaflet)
set_restoration_plot_centres <- function(input = NULL, width = NULL, length = NULL) {
# Calculate the coordinates of the rectangular polygon
x <- sf::st_coordinates(input)[1, 1]
y <- sf::st_coordinates(input)[1, 2]
# set parameters
x_min <- x - (width / 2)
x_max <- x + (width / 2)
y_min <- y - (length / 2)
y_max <- y + (length / 2)
polygon <- sf::st_polygon(list(rbind(c(x_min, y_min), c(x_min, y_max), c(x_max, y_max), c(x_max, y_min), c(x_min, y_min)))) |>
sf::st_sfc(crs = 20353)
return(polygon)
}
habitats <- c("Plateau", "Back Reef Slope", "Reef Slope", "Sheltered Reef Slope", "Inner Reef Flat", "Outer Reef Flat", "Reef Crest")
moore_geomorphic <- st_read("/Users/rof011/GBR-coral-restoration/data/Moore-Reef-20230906030257/Geomorphic-Map/geomorphic.geojson", quiet=TRUE) %>%
sf::st_transform(crs = sf::st_crs("EPSG:20353")) %>%
dplyr::group_by(class) %>%
dplyr::mutate(habitat_id = paste0(
gsub(" ", "_", class), "_",
sprintf(paste0("%0", ceiling(log10(max(1:length(class)))), "d"), 1:length(class)))) %>%
sf::st_make_valid() %>%
# mutate(habitat_id=as.factor(habitat_id)) %>%
dplyr::group_by(habitat_id, class) %>%
summarise(geometry = st_union(geometry)) %>%
mutate(area = round(st_area(geometry))) %>%
filter(class %in% habitats) %>%
drop_na(class)
moore_plots_centroids <- moore_geomorphic %>%
filter(area > set_units(10000,"m^2")) %>%
filter(class %in% c("Reef Slope", "Sheltered Reef Slope", "Back Reef Slope")) %>%
drop_na(class) %>%
st_centroid() %>%
st_cast("POINT")
moore_plot_outputs <- foreach(i=1:nrow(moore_plots_centroids), .combine="c") %do% {
tmp <- set_restoration_plot_centres(moore_plots_centroids$geometry[i], width=100, length=100)
tmp
}
generate_random_points <- function(polygon, n) {
st_sample(polygon, size = n)
}
moore_plot_outputs3 <- moore_plot_outputs |> st_as_sf() |> slice_sample(n=1)
moore_plots_centroids_small <- st_as_sf(st_centroid(moore_plot_outputs3))
#moore_plot_small <- set_restoration_plot_centres(moore_plot_outputs3, width=10, length=10) #|> st_as_sf() |> mutate(id=moore_plot_outputs3$id)
width=10
length=10
# Calculate the coordinates of the rectangular polygon
x <- sf::st_coordinates(moore_plots_centroids_small)[1] |> as.numeric()
y <- sf::st_coordinates(moore_plots_centroids_small)[2] |> as.numeric()
# set parameters
x_min <- x - (width / 2)
x_max <- x + (width / 2)
y_min <- y - (length / 2)
y_max <- y + (length / 2)
moore_plot_small <- sf::st_polygon(list(rbind(c(x_min, y_min), c(x_min, y_max), c(x_max, y_max), c(x_max, y_min), c(x_min, y_min)))) |>
sf::st_sfc(crs = 20353)
all_random_points <- st_sample(moore_plot_small, size = 5000) |> st_as_sf(crs=20353)
coral_locations <- all_random_points %>%
mutate(coral_id=seq(1,n(),1)) %>%
mutate(coral_id=as.factor(coral_id)) %>%
st_buffer(dist = 0.01)
coral_locations_yr5 <- all_random_points %>%
slice_sample(n=1000) %>%
mutate(coral_id=seq(1,n(),1)) %>%
#mutate(radius=2.5) |>
st_buffer(dist = 0.01+(0.025*4))
#moore_satellite_inset <- raster::crop(moore_satellite, st_bbox(moore_plot_outputs2 |> st_transform(3857)))
# make an internal box of 10x10 to speed things up
map <- tm_shape(moore_plot_outputs3, is.master=TRUE) +
tm_borders(col=NA) +
tm_tiles("Esri.WorldImagery", group = "<b> [Seascape]</b> satellite map", alpha = 0.6) +
tm_shape(moore_satellite_inset, raster.downsample=FALSE, alpha=0.6) +
tm_rgb(r=1, g=2, b=3) +
tm_shape(moore_geomorphic, id="area", name = "<b> [Seascape]</b> habitats") +
tm_borders(col = "black", lwd = 0.2) +
tm_fill("class", name = "Benthic habitats", id="area", alpha = 0.6) +
tm_shape(moore_plot_outputs3, name = "<b> [Restoration]</b> 1ha plot") +
tm_borders(col="darkred", lwd=1.5) +
tm_shape(moore_plot_small, name = "<b> [Restoration]</b> 100m2 plot", is.master=TRUE) +
tm_borders(col="darkred", lwd=1.5) +
tm_shape(coral_locations, id="coral_id", legend.show=FALSE, name = "<b> [Corals]</b> Year 1") +
tm_polygons(col="darkblue", alpha=0.4) +
tm_shape(coral_locations_yr5, id="coral_id", legend.show=FALSE, name = "<b> [Corals]</b> Year 5") +
tm_polygons(col="aquamarine3", alpha=0.4) +
tm_view(set.zoom.limits=c(5,100), bbox=st_bbox(moore_plot_small))
map |> tmap::tmap_leaflet() |>
leaflet::addProviderTiles('Esri.WorldImagery', group = "<b> [Seascape]</b> satellite map", options=leaflet::providerTileOptions(maxNativeZoom=18,maxZoom=100)) |>
leaflet::addProviderTiles('Esri.WorldTopoMap', group = "<b> [Seascape]</b> base map", options=leaflet::providerTileOptions(maxNativeZoom=19,maxZoom=100)) |>
leaflet::addLayersControl(position="bottomleft", overlayGroups=c(
"<b> [Seascape]</b> satellite map", "<b> [Seascape]</b> habitats", "<b> [Restoration]</b> plots",
"<b> [Restoration]</b> 1ha plot", "<b> [Restoration]</b> inset plot", "<b> [Corals]</b> Year 1", "<b> [Corals]</b> Year 5"),
options=leaflet::layersControlOptions(collapsed = FALSE)) |>
leaflet.extras::addFullscreenControl(position = "bottomleft", pseudoFullscreen = FALSE)